No evidence of neural feature-specific pre-activation during the prediction of an upcoming stimulus

A reanalysis of Demarchi et al. (Nature Communications, 2019)

Authors
Affiliation

Abdoun Oussama

Centre de Rercherches en Neurosciences de Lyon

Todorov Dimitrii

Poublan-couzardot Arnaud

Tirou Coumarane

Lutz Antoine

Vernet Marine

Quentin Romain

README

Background information

See the associated Github project home page.

The present R Quarto document calculates the correlation between decoding accuracy and sequence entropy, perform cluster-based permutation statistical tests and Bayesian statistics on accuracy and correlation time-generalization matrices, and plot all the figures and supplementary figures of the paper.

Instructions for use

  1. Obtain the classification data either from Zenodo, or by running the Python scripts on raw MEG data (see the Github project home page for details). Save it in the ./data folder.

  2. For a first quick run of the script, we recommend that you set the number of cluster permutations to a low value (e.g. 100) in Prepare / Parameters / Global. Set it to 10,000 for stable results matching the ones reported in the publication.

Prepare

Libraries

Source code
library(expm)         # exponentiate matrices
Loading required package: Matrix

Attaching package: 'expm'
The following object is masked from 'package:Matrix':

    expm
Source code
library(patchwork)    # combine plots
library(tictoc)       # time executions
library(correlation)
library(abind)        # combine multi-dimensional arrays
library(BayesFactor)
Loading required package: coda
************
Welcome to BayesFactor 0.9.12-4.6. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).

Type BFManual() to open the manual.
************
Source code
library(magrittr)     # special pipes
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.0     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand()    masks Matrix::expand()
✖ tidyr::extract()   masks magrittr::extract()
✖ dplyr::filter()    masks stats::filter()
✖ dplyr::lag()       masks stats::lag()
✖ tidyr::pack()      masks Matrix::pack()
✖ purrr::set_names() masks magrittr::set_names()
✖ tidyr::unpack()    masks Matrix::unpack()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Source code
library(reticulate)
np <- import("numpy")
os <- import("os")
mne <- import("mne")

Parameters

Global

Source code
path.root <- "./data/results"
list.subj <- list.dirs(path.root, full.names = F, recursive = F)
n.subj <- length(list.subj)

time_offset <- 670
n_permutations <- 100
save <- FALSE

Plotting

Source code
# Temporal resolution in ms (1000 / sampling frequency)
dt <- 10

# Colors
col.pal.base <- rev(pals::brewer.rdbu(10))
col.pal.diff <- rev(pals::brewer.brbg(10))
col.pal.corr <- pals::ocean.curl(100)
col.lines <- "grey40"
col.clusters <- "grey10"

# Breaks & Labels
# --- function
format_labels_acc <- function(labels) {
  labels[seq(2,length(labels),2)] <- ""
  labels <- str_replace(labels, "0.", ".")
}
# --- base
acc.breaks.base <- seq(.21,.29, .005)
acc.labels.base <- format_labels_acc(acc.breaks.base)
# --- diff
acc.breaks.diff <- seq(-.02,.02, .005)
acc.labels.diff <- format_labels_acc(acc.breaks.diff)
# --- correlations
corr.breaks <- seq(-1,1,.2)
corr.labels <- str_replace(corr.breaks, "0.", ".")

# GGPLOT theme
theme_set(theme_minimal())
theme_update(text = element_text(family = "Arial"),
             plot.tag.position = c(0, 0.95),
             plot.tag = element_text(size = 14, hjust = 0, vjust = 1, face = "bold"))

Transition matrices

From Demarchi, Sanchez, and Weisz (2019)

Source code
mat.T <- list(RD = matrix(rep(.25, 16), ncol = 4, byrow = T),
              MM = matrix(), MP = matrix(), OR = matrix())

mat.T$MM <- matrix(c(c(.25,0,.37,.38),
                     c(.38,.25,0,.37),
                     c(.37,.38,.25,0),
                     c(0,.37,.38,.25)),
                   ncol = 4, byrow = T)

mat.T$MP <- matrix(c(c(.25,0,.15,.60),
                     c(.60,.25,0,.15),
                     c(.15,.60,.25,0),
                     c(0,.15,.60,.25)),
                   ncol = 4, byrow = T)

mat.T$OR <- matrix(c(c(.25,0,0,.75),
                     c(.75,.25,0,0),
                     c(0,.75,.25,0),
                     c(0,0,.75,.25)),
                   ncol = 4, byrow = T)

covCT <- function(mat.C, mat.T) {
  covar <- 0
  for (i in 1:dim(mat.C)[1]) {
    covar <- covar + mean((mat.C[i,]-mean(mat.C[i,]))*(mat.T[i,]-mean(mat.T[i,])))
  }
  return(covar)
}

Functions

Matrix to df

Source code
matrix_to_df2 <- function(array, row.name, col.name, val.name) {
  array %>% 
    as_tibble() %>%
    rownames_to_column(row.name) %>%
    pivot_longer(-row.name, names_to = col.name, values_to = val.name) %>% 
    mutate("{col.name}" := str_remove(!!sym(col.name), "V")) %>% 
    mutate_all(as.numeric) %>% 
    return
}

# !!! USE as.data.frame.table() instead

Batch data loading

Source code
load_loop_subj <- function(root = path.root, folder, condition, diag = F, crop=T) {
  
  # Individual scores
  if (folder == "omissions") {
    array.scores <- array(dim = c(141,141,n.subj))
  } else if (folder == "sounds") {
    array.scores <- array(dim = c(33,135,n.subj))
  }
  
  i <- 1
  for (s in list.subj) {
    file <- os$path$join(root, s, folder,
                         paste0("cv_", condition, "_scores.npy"))

    # Proceed to next subject if file doesn't exist
    if (!os$path$isfile(file)) {
      cat(file, " does not exist !!!\n")
      next
    }
    
    # Load data otherwise
    tmp <- np$load(file)
    
    # Average across CV folds if present
    if (length(dim(tmp)) == 3) {
      tmp <- tmp %>% apply(c(2,3), mean) # average across CV folds
    }
      
    array.scores[,,i] <- tmp
    i <- i+1
  }
  
  # Average across participants
  t_train_offset <- case_when(folder == "omissions" ~ time_offset,
                              folder == "sounds" ~ 0)
  t_test_offset <- time_offset
  
  df.mean <- array.scores %>% 
    apply(c(1,2), mean) %>%
    matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "accuracy") %>%
    # --- convert indices to time
    mutate(t_train = dt*(t_train-1) - t_train_offset,
           t_test =  dt*(t_test-1)  - t_test_offset)
  
  # --- crop time if requested (default)
  if (crop) {
    df.mean %<>% filter(t_train>=0, (t_train-333) < dt)
  }
  
  return(list(array.scores, df.mean))
}

Correlations

Source code
## Correlation function between two vectors
correlations_vecs <- function(vec1, vec2) {
  cor(as.vector(vec1), as.vector(vec2))
}

# Bind time-generalization matrices from different entropy conditions into a single 3D array
correlations_extract3D <- function(data, idx.lines, idx.subj) {
  arr.tmp <- data[[idx.lines[1],"array"]][[1]][,,idx.subj]
  for (i in 2:length(idx.lines)) {
    arr.tmp <- abind(arr.tmp, data[[idx.lines[i],"array"]][[1]][,,idx.subj], along = 3)
  }
  return(arr.tmp)
}


# Replace NAs by 0 and extreme values by +/- .99
correlations_fixval <- function(arr.corr) {
  arr.corr[is.na(arr.corr)] <- 0
  arr.corr[arr.corr == 1] <- .99
  arr.corr[arr.corr == -1] <- -.99
  return(arr.corr)
}

Cluster permutation stats

Source code
get_signif_clust2d <- function(X, n_permutations = 2**12, n_jobs = -1) {
  # run permutaiton tests
  mne$stats$spatio_temporal_cluster_1samp_test(X,
                                               n_permutations = n_permutations,
                                               n_jobs = as.integer(n_jobs),
                                               out_type = "mask",
                                               verbose = FALSE) -> tmp
  names(tmp) <- c("t_obs","clusters","cluster_pv","H0")
  
  # Extract clusters: the final output is a 0/1 matrix, 1 indicating belonging to a significant cluster 
  # --- get ids of significant clusters
  idx.signif <- which(tmp$cluster_pv < .05) 
  # --- aggregate clusters
  signif <- array(0, dim = dim(tmp$t_obs)) 
  for (i in idx.signif) {
    signif <- signif + tmp$clusters[[i]]
  }
  
  return(1*(signif > 0))
}


# Loop through multiple conditions, stored in a dataframe where each row contains a 3D array in an "array" column
# Other columns are appended to the output as descriptors of the condition
cluster_loop_cond <- function(df, H0 = 0, n = 2**12, n_jobs = -1, crop = FALSE) {
  # Initialize significant clusters dataframe
  df.clusters <- tibble()
  # Loop through all conditions
  for (i in 1:nrow(df)) {
    
    # Extract array of subject-level time-generalized accuracy
    X <- df[[i,"array"]][[1]]
    # --- crop to stimulus N-1 if requested
    if (crop) {X <- X[,35:68,]}
    # --- reshape to how MNE expects it
    X <- X %>% np$moveaxis(as.integer(2), as.integer(0)) 
  
    df.clusters %<>% bind_rows(
      get_signif_clust2d(X-H0, n_jobs = n_jobs, n_permutations = n) %>% 
        matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "signif") %>% 
        bind_cols(select(df,-array)[i,])
    )
  }
  
  return(df.clusters)
}

Difference original-reordered

Source code
calculate_diff <- function(df, f.map, values_to) {
  df %>% 
    mutate(manip = case_match(manip, "" ~ "original", "_reord" ~ "reordered")) %>% 
    pivot_wider(names_from = manip, values_from = values_to) %>% 
    mutate(diff = f.map(original, reordered, .f=\(x,y){x-y})) %>% 
    pivot_longer(cols = c(original, reordered, diff), names_to = "manip", values_to = values_to)
}

Base plot

Source code
plot_base <- function(df.plot,
                      z = "accuracy",
                      contour = TRUE,
                      col.contour = col.clusters,
                      col.pal = rev(pals::brewer.rdbu(10)),
                      z.breaks = acc.breaks.base,
                      z.labels = acc.labels.base,
                      legend.position = "bottom") {
  
  # Onset lines
  lines.v.pos <- 330*seq(round(min(df.plot$t_test)/333), round(max(df.plot$t_test)/333), 1)
  lines.v.pos <- lines.v.pos[lines.v.pos!=0]
  lines.h.pos <- 330*seq(round(min(df.plot$t_train)/333), round(max(df.plot$t_train)/333), 1)
  lines.h.pos <- lines.h.pos[lines.h.pos!=0]

  # Plot
  g <- df.plot %>% 
    ggplot(aes(x = t_test, y = t_train)) +
    # --- heatmap
    geom_tile(aes_string(fill = z)) +
    # --- sounds onsets
    geom_vline(xintercept = lines.v.pos, linetype = 2, color = col.lines) +
    geom_hline(yintercept = lines.h.pos, linetype = 2, color = col.lines) +
    geom_vline(xintercept = 0, color = col.lines) +
    geom_hline(yintercept = 0, color = col.lines) +
    
    scale_x_continuous(breaks = scales::pretty_breaks(12),
                       expand = c(0,0)) +
    scale_y_continuous(breaks = scales::pretty_breaks(4),
                       expand = c(0,0)) +
    scale_fill_stepsn(colors = col.pal,
                      breaks = z.breaks,
                      labels = z.labels,
                      limits = c(min(z.breaks), max(z.breaks)),
                      oob = scales::squish) +
    
    labs(x = "Test time (ms)", y = "Train time (ms)") +
  
    coord_equal(clip = "off") +
    
    theme(# Text size
          plot.title = element_text(size = 12, hjust = 0.5, face = "bold"),
          plot.subtitle = element_text(size = 9, hjust = 0.5),
          axis.title = element_text(size = 9),
          legend.title = element_text(size = 9),
          strip.text = element_text(size = 9),
          axis.text = element_text(size = 5),
          legend.text = element_text(size = 5.5),
          # Time ticks
          axis.ticks = element_line(color = "black"),
          axis.ticks.length = unit(0.1, "cm"),
          # Legend
          legend.position = legend.position,
          legend.box.margin = margin(0,0,0,0),
          legend.box.spacing = unit(0, "pt"),
          # no grid
          panel.grid = element_blank(),
          # don't clip facet titles
          strip.clip = "off")
  
  if (legend.position %in% c("bottom","top")) {
    g <- g + 
      guides(fill = guide_colorbar(barwidth = 0.5*length(z.breaks-1), barheight = 0.5,
                                   ticks = F, title.vjust = 1))
  } else {
    g <- g + 
      guides(fill = guide_colorbar(barheight = 0.5*length(z.breaks-1), barwidth = 0.5,
                                   ticks = F, title.vjust = 1))
  }
  
  # --- significant clusters
  if (contour) {
    g <- g +
      geom_contour(aes(z = signif), size = 0.1, color = col.contour)
  }
  
  return(g)
}

Replication & Reordering (empirical null)

Load data

Sounds

Source code
tic()

# Load subject-level data (used for clustering) and derive grand average (used for visualization)
df.stim.mean <- tibble()
df.stim.subj <- tibble()
for (manip in c("","_reord")) { 
  
  for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
    tmp <- load_loop_subj(folder = "sounds",
                          condition = paste0(direction, manip))
    
    df.stim.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
    df.stim.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  }
}
Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.
Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(row.name)

  # Now:
  data %>% select(all_of(row.name))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
Source code
df.stim.mean %<>% mutate(stim = "sounds")
df.stim.subj %<>% mutate(stim = "sounds") %>%
  # --- center decoding accuracy on 0
  mutate(array = map(array, .f=\(x){x-0.25}))

toc() # ~1s
1.253 sec elapsed

Omissions

Source code
tic()

df.omit.mean <- tibble()
df.omit.subj <- tibble()
for (manip in c("","_reord")) { 
  
  for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
      
    tmp <- load_loop_subj(folder = "omissions",
                          condition = paste0(direction, manip))
    
    df.omit.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
    df.omit.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  }
}
Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.
Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(row.name)

  # Now:
  data %>% select(all_of(row.name))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
./data/results/11920812IONP/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19750430PNRK/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19800616MRGU/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19810726GDZN/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19810918SLBR/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19821223KRHR/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19830114RFTM/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19851130EIFI/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19861231ANSR/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19871229CRBC/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19880328AGSG/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19891222GBHL/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19900606KTAD/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19901026KRKE/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19910823SSLD/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19920804CRLE/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19921111BRHC/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930118IMSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930423EEHB/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930524GNTA/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930621ATLI/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930630MNSU/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19940403HIEC/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19940601IGSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19940930SZDT/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19950212BTKC/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19950326IIQI/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19950604LLZM/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19950905MCSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19960304SBPE/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19960418GBSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19960708HLHY/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19961118BRSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
Source code
df.omit.mean %<>% mutate(stim = "omissions") %>% filter(abs(t_test) <= 670)
df.omit.subj %<>% mutate(stim = "omissions") %>%
  # --- center decoding accuracy on 0
  mutate(array = map(array, .f=\(x){x-0.25})) %>% 
  # --- crop subject-level data for cluster analysis
  mutate(array = map(array, .f=\(x){x[71:104,4:138,]}))

toc() # ~ 2s
3.59 sec elapsed

Correlations with entropy

Source code
## Prepare input data
df.tmp <- bind_rows(df.stim.subj, df.omit.subj) %>%
  calculate_diff(f.map = map2, values_to = "array") %>%
  mutate(regularity = str_remove(direction, "rd_to_"),
         regularity = factor(regularity, levels = c("rd","mm","mp","or"))) %>% 
  select(-direction) %>% 
  arrange(as.numeric(regularity))
df.tmp
# A tibble: 24 × 4
   stim      manip     array                 regularity
   <chr>     <chr>     <list>                <fct>     
 1 sounds    original  <dbl [33 × 135 × 33]> rd        
 2 sounds    reordered <dbl [33 × 135 × 33]> rd        
 3 sounds    diff      <dbl [33 × 135 × 33]> rd        
 4 omissions original  <dbl [34 × 135 × 33]> rd        
 5 omissions reordered <dbl [34 × 135 × 33]> rd        
 6 omissions diff      <dbl [34 × 135 × 33]> rd        
 7 sounds    original  <dbl [33 × 135 × 33]> mm        
 8 sounds    reordered <dbl [33 × 135 × 33]> mm        
 9 sounds    diff      <dbl [33 × 135 × 33]> mm        
10 omissions original  <dbl [34 × 135 × 33]> mm        
# ℹ 14 more rows
Source code
tic()
df.corr.mean <- tibble()
df.corr.subj <- tibble() 

## Loop through condition "manip" and stimulus type "stim
for (mm in c("original","reordered","diff")) { 

  for (ss in c("sounds", "omissions")) { #
    cat("\n\nProcessing ", str_to_upper(mm), " / ", str_to_upper(ss), "...\n", sep="")
    
    # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
    for (idx.subj in 1:n.subj) {
      
      # --- bind time-generalization matrices for the 3 non-random conditions into a single 3D array
      # (add accuracy data from the random condition if we are processing sounds)
      if (ss == "sounds") {
        idx.lines = 1:4
      } else if (ss == "omissions") {
        idx.lines = 2:4
      }
      arr.tmp <- correlations_extract3D(data = filter(df.tmp, manip == mm, stim == ss),
                                        idx.lines = idx.lines, idx.subj)
      
      # --- calculate correlation with entropy level coded as 0,1,2(,3), obtaining a 2D array
      if (ss == "sounds") {
        arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:3)
      } else if (ss == "omissions") {
        arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 1:3)
      }
      
      # --- append to already processed subjects in a 3D array
      if (idx.subj == 1) {
        arr.corr <- arr.corr.tmp
      } else {
        arr.corr %<>% abind(arr.corr.tmp, along = 3)
      }
    }
    
    # Postprocessing
    arr.corr %<>% 
      correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
      atanh # apply Fisher's transformation
    
    # --- store participant-level results in a dataframe
    df.corr.subj %<>% bind_rows(tibble(manip = mm, stim = ss, array = list(arr.corr)))
    
    # --- calculate groupe average and store in a dataframe
    df.corr.mean %<>% bind_rows(
      arr.corr %>% 
        apply(c(1,2), mean) %>% # average across participants
        tanh %>% # convert back to correlation coefficients  
        matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>% 
        mutate(manip = mm, stim = ss)
    )
  }
}


Processing ORIGINAL / SOUNDS...


Processing ORIGINAL / OMISSIONS...


Processing REORDERED / SOUNDS...


Processing REORDERED / OMISSIONS...


Processing DIFF / SOUNDS...


Processing DIFF / OMISSIONS...
Source code
toc() # ~20 s
33.512 sec elapsed

Clusters stats

Accuracy

Source code
tic()

df.tmp.subj <- bind_rows(df.stim.subj, df.omit.subj)
df.tmp.mean <- bind_rows(df.stim.mean, df.omit.mean)
  
# Cluster analysis
df.clusters.0 <- df.tmp.subj %>% 
  calculate_diff(f.map = map2, values_to = "array") %>% 
  cluster_loop_cond(n = n_permutations) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))
Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(values_to)

  # Now:
  data %>% select(all_of(values_to))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.
Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(row.name)

  # Now:
  data %>% select(all_of(row.name))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
Source code
# Quick check of clusters with a mask plot
df.clusters.0 %>% 
  mutate(direction = (direction %>% str_replace("rd_to_", "") %>% str_to_upper)) %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid(direction ~ manip) +
  geom_tile() +
  coord_equal()

Source code
# Aggregate all data and clean
df.plot.0 <- df.tmp.mean %>%
  calculate_diff(f.map = map2_dbl, values_to = "accuracy") %>% 
  # --- append cluster significance
  left_join(df.clusters.0) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
                                diff = "ORIGINAL – REORDERED"))) %>% 
  mutate(direction = toupper(sub("rd_to_", "", direction)),
         direction = factor(direction, levels = rev(c("RD","MM","MP","OR"))))
Joining with `by = join_by(t_train, t_test, direction, stim, manip)`
Source code
toc() # ~3.5 minutes
23.236 sec elapsed

Correlations

Source code
tic()

# Cluster analysis
df.clusters.corr <- df.corr.subj %>% 
  cluster_loop_cond(n = n_permutations)

# Quick check of clusters with a mask plot
df.clusters.corr %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid( ~ manip) +
  geom_tile() +
  coord_equal()

Source code
# --- append cluster significance
df.plot.0.corr <- df.corr.mean %>% 
  left_join(df.clusters.corr) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
                                diff = "ORIGINAL – REORDERED"))) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))
Joining with `by = join_by(t_train, t_test, manip, stim)`
Source code
df.plot.0.corr %>% summary
    t_train          t_test           r                             manip     
 Min.   :  0.0   Min.   :-670   Min.   :-0.91933   ORIGINAL            :9045  
 1st Qu.: 80.0   1st Qu.:-340   1st Qu.:-0.21316   REORDERED           :9045  
 Median :160.0   Median :   0   Median :-0.01320   ORIGINAL – REORDERED:9045  
 Mean   :162.5   Mean   :   0   Mean   :-0.01764                              
 3rd Qu.:250.0   3rd Qu.: 340   3rd Qu.: 0.18399                              
 Max.   :330.0   Max.   : 670   Max.   : 0.86931                              
     stim               signif      
 Length:27135       Min.   :0.0000  
 Class :character   1st Qu.:0.0000  
 Mode  :character   Median :0.0000  
                    Mean   :0.1094  
                    3rd Qu.:0.0000  
                    Max.   :1.0000  
Source code
toc() # ~ 1 min
3.179 sec elapsed

Theoretical null

Load confusion matrices & calculate

Source code
tic()

null_accuracy <- function(mat.confusion, mat.transition) {
  # Calculate accuracy under the null hypothesis from the transition and confusion matrices mat.T and mat.C
  # The 3/4 factor is used to compensate the unbiased estimator used in cov() by default (see ?cov)
  return(0.25 + sum(diag(3/4*cov(t(mat.transition), t(mat.confusion)))))
}

df.maths.subj <- tibble()
df.maths.mean <- tibble()
for (entropy in c("RD","MM","MP","OR")) {

  array.scores <- array(0.25, dim = c(33,136,n.subj))
  
  i <- 1
  for (s in list.subj) {
    
    # --- load the confusion matrix
    tmp.file <- os$path$join(path.root, s, "rd_to_rd_confmats.npz") %>%
      np$load()
    
    # --- crop around S0 (decoded stimulus)
    mat.C <- tmp.file$f[["arr_0"]][,71:103,71:104,,] %>%
      apply(c(2,3,4,5), mean)
    
    # NOTE: because the transition matrix is the same for all subjects, theoretical accuracy can be calculated from the group average of the confusion matrix (faster)
    # --- calculate theoretical accuracy
    #   --- for S-2
    array.scores[1:33,1:34,i] <- mat.C %>%
      apply(c(1,2), null_accuracy, mat.transition = mat.T[[entropy]]%^%2)
    #   --- for S-1
    array.scores[1:33,35:68,i] <- mat.C %>%
      apply(c(1,2), null_accuracy, mat.transition = mat.T[[entropy]]%^%1)
    #   --- for S
    array.scores[1:33,69:102,i] <- mat.C %>%
      apply(c(1,2), \(x)mean(diag(x)))
    #   --- for S+1
    array.scores[1:33,103:136,i] <- mat.C %>%
      apply(c(1,2), null_accuracy, mat.transition = t(mat.T[[entropy]]))
    i <- i+1
  }
  
  array.scores <- array.scores[,2:136,]
  
  # Append to other entropy levels in dataframe
  # --- subject-level data
  df.maths.subj %<>% bind_rows(tibble(array = list(array.scores), entropy = entropy))
  # --- group-averaged data
  df.maths.mean %<>% bind_rows(array.scores %>% 
      apply(c(1,2), mean) %>%
      matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "accuracy") %>% 
    mutate(entropy = entropy)
  )
}

toc() # ~ 1min
67.387 sec elapsed

Append original data & calculate diff

Source code
df.maths.all <- 
# --- aggregate original, theoretical and difference between the 2
  bind_rows(df.stim.subj %>%
              # --- select original analysis
              filter(manip == "") %>% mutate(manip = "original") %>% 
              # --- crop to the same time-window than for theoretical analysis
              mutate(entropy = toupper(str_replace(direction,"rd_to_",""))),
            df.maths.subj %>%
              mutate(manip = "maths") %>% 
              mutate(array = map(array, .f=\(x){x-0.25}))) %>% 
  select(manip, entropy, array) %>% 
  pivot_wider(names_from = manip, values_from = "array") %>% 
  mutate(diff = map2(original, maths, .f=\(x,y){x-y})) %>% 
  pivot_longer(cols = c(original, maths, diff), names_to = "manip", values_to = "array")

Correlations with entropy

Source code
tic()

df.corr.maths.mean <- tibble()
df.corr.maths.subj <- tibble() 

for (mm in c("original","maths","diff")) { 
  # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
  for (idx.subj in 1:n.subj) {
    
    # --- bind time-generalization matrices of the 4 conditions into a single 3D array
    arr.tmp <- correlations_extract3D(data = df.maths.all %>% filter(manip == mm),
                                      idx.lines = 1:4, idx.subj)
    
    # --- calculate correlation with entropy level coded as 0,1,2,3, obtaining a 2D array
    arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:3)
    
    # --- append to already processed subjects in a 3D array
    if (idx.subj == 1) {
      arr.corr <- arr.corr.tmp
    } else {
      arr.corr %<>% abind(arr.corr.tmp, along = 3)
    }
  }
  
  # Postprocessing
  arr.corr %<>% 
    correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
    atanh # apply Fisher's transformation
    
  # --- store participant-level results in a dataframe of arrays
  df.corr.maths.subj %<>% bind_rows(tibble(manip = mm, array = list(arr.corr)))
  
  # --- calculate group average and store in a dataframe
  df.corr.maths.mean %<>% bind_rows(
    arr.corr %>% 
      apply(c(1,2), mean) %>% # average across participants
      tanh %>% # convert back to correlation coefficients  
      matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>% 
      mutate(manip = mm)
  )
}

toc() # ~ 20s
16.723 sec elapsed

Clusters stats

Accuracy

Source code
tic()

df.clusters.maths <- df.maths.all %>% 
  cluster_loop_cond(n = n_permutations) %>%
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))

toc() # ~ 2 min
2.065 sec elapsed
Source code
# Quick check of clusters with a mask plot
df.clusters.maths %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid(entropy ~ manip) +
  geom_tile() +
  coord_equal()

Correlations

Source code
tic()
df.clusters.corr.maths <- df.corr.maths.subj %>% 
  cluster_loop_cond(n = n_permutations)
toc() # ~ 20s
0.566 sec elapsed
Source code
# Quick check of clusters with a mask plot
df.clusters.corr.maths %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid( ~ manip) +
  geom_tile() +
  coord_equal()

Source code
# --- append cluster significance
df.plot.maths.corr <- df.corr.maths.mean %>% 
  left_join(df.clusters.corr.maths) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", maths = "THEORETICAL NULL",
                                diff = "ORIGINAL – THEORETICAL NULL"))) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))
Joining with `by = join_by(t_train, t_test, manip)`
Source code
df.plot.maths.corr %>% summary
    t_train        t_test           r           
 Min.   :  0   Min.   :-670   Min.   :-0.98948  
 1st Qu.: 80   1st Qu.:-340   1st Qu.:-0.23120  
 Median :160   Median :   0   Median : 0.00000  
 Mean   :160   Mean   :   0   Mean   :-0.03176  
 3rd Qu.:240   3rd Qu.: 340   3rd Qu.: 0.16949  
 Max.   :320   Max.   : 670   Max.   : 0.95706  
                         manip          signif      
 ORIGINAL                   :4455   Min.   :0.0000  
 THEORETICAL NULL           :4455   1st Qu.:0.0000  
 ORIGINAL – THEORETICAL NULL:4455   Median :0.0000  
                                    Mean   :0.1631  
                                    3rd Qu.:0.0000  
                                    Max.   :1.0000  

Prediction of the most likely (“simple prediction”)

Load data

Source code
tic()

# Load subject-level data and derive grand average
df.stim.ml.mean <- tibble()
df.stim.ml.subj <- tibble()
for (manip in c("","_reord")) { 
  
  for (direction in c("rd_to_mm", "rd_to_mp", "rd_to_or")) {
    
    tmp <- load_loop_subj(folder = "sounds",
                          condition = paste0(direction, manip, "_sp"))
    
    df.stim.ml.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
    df.stim.ml.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  }
}

df.stim.ml.mean %<>% mutate(stim = "sounds")
df.stim.ml.subj %<>% mutate(stim = "sounds") %>%
  # --- center decoding accuracy on 0
  mutate(array = map(array, .f=\(x){x-0.25}))

toc() # ~ 1s
0.891 sec elapsed
Source code
df.stim.ml.mean %>% mutate_if(is.character, as.factor) %>% summary
    t_train        t_test        accuracy         manip          direction   
 Min.   :  0   Min.   :-670   Min.   :0.2040         :13365   rd_to_mm:8910  
 1st Qu.: 80   1st Qu.:-340   1st Qu.:0.2465   _reord:13365   rd_to_mp:8910  
 Median :160   Median :   0   Median :0.2498                  rd_to_or:8910  
 Mean   :160   Mean   :   0   Mean   :0.2500                                 
 3rd Qu.:240   3rd Qu.: 340   3rd Qu.:0.2532                                 
 Max.   :320   Max.   : 670   Max.   :0.3095                                 
     stim      
 sounds:26730  
               
               
               
               
               

Correlations with entropy

Source code
tic()

## Prepare input data
df.tmp <- df.stim.ml.subj %>% 
  calculate_diff(f.map = map2, values_to = "array") %>%
  mutate(entropy = str_remove(direction, "rd_to_"),
         entropy = factor(entropy, levels = c("rd","mm","mp","or"))) %>% 
  select(-direction) %>% 
  arrange(as.numeric(entropy))


## Calculate correlations
df.corr.ml.mean <- tibble()
df.corr.ml.subj <- tibble() 

# --- Loop through condition "manip" and stimulus type "stim"
for (mm in c("original","reordered","diff")) { 

  for (ss in c("sounds")) {
    cat("\n\nProcessing ", str_to_upper(mm), " / ", str_to_upper(ss), "...\n", sep="")

    
    # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
    for (idx.subj in 1:n.subj) {
      
      # --- bind time-generalization matrices for the 3 non-random conditions into a single 3D array
      arr.tmp <- correlations_extract3D(data = filter(df.tmp, manip == mm, stim == ss),
                                        idx.lines = 1:3, idx.subj)
      
      # --- calculate correlation with entropy level coded as 0,1,2, obtaining a 2D array
      arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:2)
      
      # --- append to already processed subjects in a 3D array
      if (idx.subj == 1) {
        arr.corr <- arr.corr.tmp
      } else {
        arr.corr %<>% abind(arr.corr.tmp, along = 3)
      }
    }
    
    # Post-processing
    arr.corr <- arr.corr %>% 
      correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
      atanh # apply Fisher's transformation
    
    # --- store participant-level results in a dataframe
    df.corr.ml.subj %<>% bind_rows(tibble(manip = mm, stim = ss, array = list(arr.corr)))
    
    # --- calculate groupe average and store in a dataframe
    df.corr.ml.mean %<>% bind_rows(
      arr.corr %>% 
        apply(c(1,2), mean) %>% # average across participants
        tanh %>% # convert back to correlation coefficients  
        matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>% 
        mutate(manip = mm, stim = ss)
    )
  }
}


Processing ORIGINAL / SOUNDS...


Processing REORDERED / SOUNDS...


Processing DIFF / SOUNDS...
Source code
toc() # ~ 15s
11.529 sec elapsed

Cluster stats

Accuracy

Source code
tic()

## Cluster analysis with subject-level data
# --- calculate original-reordered difference
tmp.diff <- df.stim.ml.subj %>% 
  calculate_diff(f.map = map2, values_to = "array")

# --- combine & run cluster analysis
df.clusters.ml <- tmp.diff %>% #bind_rows(tmp.diff, tmp.ave) %>% 
  cluster_loop_cond(n = n_permutations) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))

# Quick check of clusters with a mask plot
df.clusters.ml %>% 
  mutate(direction = (direction %>% str_replace("rd_to_", "") %>% str_to_upper)) %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid(direction ~ manip) +
  geom_tile() +
  coord_equal()

Source code
## Grand average data
# --- calculate original-reordered difference
tmp.diff <- df.stim.ml.mean %>%
  calculate_diff(f.map = map2_dbl, values_to = "accuracy")
# --- calculate average of difference across entropy levels
tmp.ave <- tmp.diff %>% 
  filter(manip == "diff") %>% 
  group_by(t_train, t_test, manip, stim) %>% 
  summarise(accuracy = mean(accuracy, na.rm=T)) %>% 
  mutate(direction = "average of\n(MM,MP,OR)")
`summarise()` has grouped output by 't_train', 't_test', 'manip'. You can
override using the `.groups` argument.
Source code
# --- combine
df.plot.ml <- tmp.diff #bind_rows(tmp.diff, tmp.ave)

##  Aggregate all data and clean
df.plot.ml %<>%
  # --- append cluster significance
  left_join(df.clusters.ml) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
                                diff = "ORIGINAL – REORDERED"))) %>%
  mutate(direction = toupper(sub("rd_to_", "", direction))) %>% 
  mutate(direction = factor(direction, levels = rev(c("AVERAGE OF\n(MM,MP,OR)","MM","MP","OR"))))
Joining with `by = join_by(t_train, t_test, direction, stim, manip)`
Source code
toc() # ~1.5min
2.544 sec elapsed

Correlations

Source code
tic()

df.clusters.corr.ml <- df.corr.ml.subj %>% 
  cluster_loop_cond(n = n_permutations)

toc() # ~ 30s
0.583 sec elapsed
Source code
# Quick check of clusters with a mask plot
df.clusters.corr.ml %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid( ~ manip) +
  geom_tile() +
  coord_equal()

Source code
# --- append cluster significance
df.plot.ml.corr <- df.corr.ml.mean %>% 
  filter(between(t_test,35,68)) %>% 
  left_join(df.clusters.corr.ml) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
                                diff = "ORIGINAL – REORDERED"))) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))
Joining with `by = join_by(t_train, t_test, manip, stim)`
Source code
df.plot.ml.corr %>% summary
    t_train        t_test           r                             manip     
 Min.   :  0   Min.   :-330   Min.   :-0.75166   ORIGINAL            :1122  
 1st Qu.: 80   1st Qu.:-250   1st Qu.:-0.24421   REORDERED           :1122  
 Median :160   Median :-165   Median :-0.05519   ORIGINAL – REORDERED:1122  
 Mean   :160   Mean   :-165   Mean   :-0.05098                              
 3rd Qu.:240   3rd Qu.: -80   3rd Qu.: 0.13428                              
 Max.   :320   Max.   :   0   Max.   : 0.77497                              
     stim               signif         
 Length:3366        Min.   :0.0000000  
 Class :character   1st Qu.:0.0000000  
 Mode  :character   Median :0.0000000  
                    Mean   :0.0008913  
                    3rd Qu.:0.0000000  
                    Max.   :1.0000000  

Figures

Methodological plots

Transition matrices (fig. 1a)

Source code
# Convert & combine transition matrices into a dataframe
df.plot <- tibble()
for (condition in c("RD","MM","MP","OR")) {
  df.plot %<>% bind_rows(
    mat.T[[condition]] %>% 
      matrix_to_df2(row.name = "from", col.name = "to", val.name = "p") %>% 
      mutate(condition = condition)
  )  
}
df.plot %<>%
  mutate(condition = factor(condition, names(mat.T)))

df.plot.pitch <- bind_rows(tibble(x = 0.1, y = 1:4, pitch = 1:4),
                           tibble(y = 0.1, x = 1:4, pitch = 1:4))
# Plot
df.plot %>% 
  ggplot(aes(x = to, y = from)) +
  facet_wrap(~ condition, nrow = 1, strip.position = "bottom") +
  # --- heatmap
  geom_tile(aes(fill = p), color = "black") +
  # --- display percentages
  geom_text(data = . %>% filter(p<=.25),
            aes(label = paste0(100*p,"%")),
            color = "black", size = 2) +
  geom_text(data = . %>% filter(p>.25),
            aes(label = paste0(100*p,"%")),
            color = "white", size = 2) +
  
  # --- display color-coded pitch in axes
  geom_point(data = df.plot.pitch,
             aes(x, y, color = pitch), size = 4, show.legend = FALSE) +
  geom_text(data = df.plot.pitch,
             aes(x, y, label = LETTERS[pitch]), color = "white", size = 2, show.legend = FALSE) +
  
  scale_color_viridis_c(end = 0.8) +
  scale_x_continuous(name = "to\n", position = "top") +
  scale_y_reverse(name = "from\n") +
  scale_fill_gradient(name = "transition probability",
                      low = "white", high = "black",
                      labels = ~str_replace(., "^0.","."),
                      limits = c(0,1)) +
  guides(fill = guide_colorbar(barheight = 0.5,
                               ticks = F, title.vjust = 1)) +
  coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
              expand = F, clip = "off") +
  theme(legend.position = "bottom",
        legend.box.margin = margin(t = -20),
        # text size
        axis.title = element_text(size = 9),
        legend.title = element_text(size = 9),
        strip.text = element_text(size = 12),
        # remove axis text
        axis.text = element_blank(),
        # Positions & margins
        axis.title.x = element_text(hjust = 0.08),
        panel.spacing.x = unit(0.8,"cm")) -> g.transitions

g.transitions

Reordering (fig. 1b)

Source code
# Generate Markov sequences
n.samples <- 12
stims <- list(RD = c(), OR = c())
set.seed(128764) #20476, 99713
for (entropy in c("RD","OR")) {
  stims[[entropy]] <- vector("numeric", n.samples)
  stims[[entropy]][1] <- sample(1:4, 1)
  for (i in 2:n.samples) {
    stims[[entropy]][i] <- sample(1:4, 1, prob = mat.T[[entropy]][stims[[entropy]][i-1],])
  }
}

# Add global sequence indices, class-specific indices and offset according to pitch
df.plot <- stims %>% as.tibble() %>% 
  mutate(idx = as.double(1:n.samples)) %>% 
  pivot_longer(-idx, names_to = "entropy", values_to = "class") %>% 
  mutate(entropy = factor(entropy, levels = c("OR","RD"))) %>%
  group_by(entropy,class) %>% mutate(rank = rank(idx)) %>% 
  mutate(y = as.double(entropy) + class/16-10/64)

pitch.y <- sort(unique(df.plot$y))

# Plot
df.plot %>% 
  ggplot(aes(x = idx, y = y, color = class)) +
  
  # --- pitch lines
  geom_hline(yintercept = pitch.y, color = "grey90") +
  # --- pitch labels
  annotate(geom = "text", label = "pitch",
           x = 12.5, y = c(mean(head(pitch.y,4)), mean(tail(pitch.y,4))),
           hjust = 0.5, vjust = 2.5, size = 2.5, color = "grey50", angle = 90) +
  
  annotate(geom = "text", label = "low",
           x = 12.5, y = c(min(head(pitch.y,4)), min(tail(pitch.y,4))),
           hjust = 0.74, vjust = 1.2, size = 2, color = "grey50", angle = 90) +
  
  annotate(geom = "text", label = "high",
           x = 12.5, y = c(max(head(pitch.y,4)), max(tail(pitch.y,4))),
           hjust = 0.25, vjust = 1.2, size = 2, color = "grey50", angle = 90) +
  
  # --- reordering trajectories
  ggh4x::geom_pointpath(data = df.plot %>% arrange(desc(entropy)),
               aes(x = idx, y = y, group = interaction(class,rank)),
               mult = 0.4, color = "black", linewidth = 0.3,
               arrow = arrow(type = "closed", angle = 25, length=unit(.15, 'cm'))) +
  # --- stimuli
  geom_point(size = 4, show.legend = F) +
  # --- stimuli ranks
  geom_text(aes(label = rank), size = 2, color = "white") +
  
  scale_x_continuous(name = "trial number", breaks = 1:n.samples) +
  scale_y_continuous(breaks = c(1,2), labels = c("OR","RD"), position = "left") +
  scale_color_viridis_c(end = 0.8) +
  coord_cartesian(clip = "off") +
  
  theme(axis.title.y = element_blank(),
        axis.text.y = element_text(size = 10, color = "black", margin = margin(r=-2)),
        axis.title.x = element_text(size = 9, color = "black"),
        panel.grid = element_blank(),
        panel.grid.major.x = element_line(color = "grey80", linewidth = 0.2)) -> g.reordering

g.reordering

Prediction types (fig. 2a)

Presented vs. most likely

Source code
# Generate Markov sequences
n.samples <- 6
stims <- mat.T
set.seed(12871) 
for (entropy in names(stims)) {
  stims[[entropy]] <- vector("numeric", n.samples)
  stims[[entropy]][1] <- 4#sample(1:4, 1)
  for (i in 2:n.samples) {
    stims[[entropy]][i] <- sample(1:4, 1, prob = mat.T[[entropy]][stims[[entropy]][i-1],])
  }
}

# Add sequence index & predictions
df.plot <- stims %>% as.tibble() %>% 
  mutate(idx = as.double(1:n.samples)) %>% 
  pivot_longer(-idx, names_to = "entropy", values_to = "class") %>% 
  mutate(entropy = factor(entropy, levels = rev(names(stims)))) %>% 
  group_by(entropy) %>% mutate(actual = lead(class),
                               mostlikely = ifelse(class==1, 4, class-1)) %>% 
  pivot_longer(c(actual,mostlikely), names_to = "status", values_to = "class.pred") %>% 
  filter(entropy != "RD")

# Plot
df.plot %>% 
  ggplot(aes(x = idx, y = class)) +
  facet_wrap(~ entropy, ncol = 1, strip.position = "right") +
  
  geom_segment(data = df.plot %>% filter(status == "mostlikely", idx != n.samples),
               aes(x = idx, xend = idx+1, y = class, yend = class.pred, linetype = as.factor(sign(class.pred-class))),
               #arrow = arrow(type = "closed", length = unit(6,"pt")),
               color = "grey50", show.legend = F) +
  
  geom_point(aes(fill = class), size = 3, shape = 21, color = "white", show.legend = T) +
  
  geom_point(data = df.plot %>% filter(status == "mostlikely", idx != n.samples),
             aes(x = idx+1, y = class.pred, color = class.pred),
             size = 2.5, shape = 8) +
  
  scale_x_continuous(name = "trial number", breaks = 1:n.samples) +
  scale_y_continuous(name = "") +
  scale_color_viridis_c(end = 0.8, guide = "legend", name = "", labels = c("","","","most likely")) +
  scale_fill_viridis_c(end = 0.8, guide = "legend", name = "", labels = c("","","","actually presented")) +
  
  guides(linetype = F) +
  coord_cartesian(clip = "off") +
  expand_limits(x = c(0.9,6.1), y = c(0.3,4.7)) +

  guides(#color = guide_legend(order = 1, override.aes = list(color = "grey50")),
         color = guide_legend(order = 1),
         fill = guide_legend(order = 2)) +
  theme(plot.caption = element_text(size = 6),
        axis.text.y = element_blank(),
        axis.title.x = element_text(size = 9, color = "black"),
        strip.text = element_text(size = 12),
        panel.spacing = unit(0.2,"cm"),
        panel.grid = element_blank(),
        panel.border = element_rect(fill=NA),
        legend.position = "top",
        legend.box = "vertical",
        legend.box.just = "left",
        legend.margin = margin(b = -10),
        legend.text = element_text(margin = margin(r = -12, l=0,  unit = "pt"))
        ) -> g.mostlikely

g.mostlikely

Figure 1 & Supplementary Figure 2

Initialize

Source code
stims <- list(RD = c(), OR = c())
g.base <- list(sound = c(), omission = c())
g.diff <- list(sound = c(), omission = c())

# Create ad hoc dataframes for onset' labels
df.labels.onset <- list(
  sounds = tibble(label = "0 = sound onset",
                manip = factor(c("REPLICATION\n", "EMPIRICAL NULL\n",
                                 "REPLICATION – EMPIRICAL NULL\n")),
                direction = factor("OR", levels = names(stims))),
  omissions = tibble(label = "0 = omission onset",
                manip = factor(c("REPLICATION\n", "EMPIRICAL NULL\n",
                                 "REPLICATION – EMPIRICAL NULL\n")),
                direction = factor("OR", levels = names(stims)))
)
  
# Function to add onset' labels to a plot
add_onset <- function(g, df.labels) {
  g <- g + geom_text(data = df.labels %>% filter(manip %in% unique(g$data$manip)),
                     aes(label = label, x = 0, y = 350),
                     hjust = 0, vjust = -0.5, size = 3, fontface = "italic")
  return(g)
}

Original & empirical null: accuracy plots

Source code
for (s in c("sounds","omissions")) { #
  # ORIGINAL & REORDERED
  df.plot.0 %>% 
    filter(stim == s, manip != "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip, 
                                 ORIGINAL = "REPLICATION\n", 
                                 REORDERED = "EMPIRICAL NULL\n")) %>%
    plot_base(col.pal = col.pal.base,
              z.breaks = acc.breaks.base,
              z.labels = acc.labels.base) +
    facet_grid(direction ~ manip) -> g
  
  g.base[["acc"]][[s]] <- g %>% add_onset(df.labels.onset[[s]])
  
  # DIFFERENCE
  g.diff[["acc"]][[s]] <- df.plot.0 %>%
    filter(stim == s, manip == "ORIGINAL – REORDERED") %>%
    mutate(manip = recode_factor(manip,
                                 `ORIGINAL – REORDERED` = "REPLICATION – EMPIRICAL NULL\n")) %>%
    plot_base(col.pal = col.pal.diff,
              z.breaks = acc.breaks.diff,
              z.labels = acc.labels.diff) +
    facet_grid(direction ~ manip) +
    theme(axis.title.y = element_blank()) -> g
  
  g.diff[["acc"]][[s]] <- g %>% add_onset(df.labels.onset[[s]])
}

# Assemble Figure 1c
g.fig1.1.acc <- (
  g.base[["acc"]][["sounds"]] + 
    (g.diff[["acc"]][["sounds"]] + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1)) &
    theme(strip.text.x = element_text(face = "bold"))
)

g.fig1.1.acc

Source code
# Assemble Supplementary Figure 1a
g.supfig2.acc <- (
  g.base[["acc"]][["omissions"]] + 
    (g.diff[["acc"]][["omissions"]] + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1)) &
    theme(strip.text.x = element_text(face = "bold"))
)

g.supfig2.acc

Original & empirical null: correlations plots

Source code
for (s in c("sounds", "omissions")) {
  # ORIGINAL & REORDERED
  df.plot.0.corr %>% 
    filter(stim == s, manip != "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip, 
                                 ORIGINAL = "REPLICATION\n", #  (original dataset)\n
                                 REORDERED = "EMPIRICAL NULL\n")) %>%  # (reordered RD)\n
    mutate(direction = "") %>% 
    plot_base(col.pal = col.pal.corr,
              z = "r",
              z.breaks = corr.breaks,
              z.labels = corr.labels) +
    labs(fill = "Correlation") +
    facet_grid(direction ~ manip) -> g
  
  g.base[["corr"]][[s]] <- g #%>% add_onset(df.labels.onset[[s]] %>% mutate(direction = ""))
  
  # DIFFERENCE
  df.plot.0.corr %>%
    filter(stim == s, manip == "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip,
                                 `ORIGINAL – REORDERED` = "REPLICATION – EMPIRICAL NULL\n")) %>% 
    mutate(direction = "") %>% 
    plot_base(col.pal = col.pal.corr,
              z = "r",
              z.breaks = corr.breaks,
              z.labels = corr.labels) +
    facet_grid(. ~ manip) +
    labs(fill = "Correlation") +
    theme(axis.title.y = element_blank()) -> g
  
  g.diff[["corr"]][[s]] <- g #%>% add_onset(df.labels.onset[[s]] %>% mutate(direction = ""))
}

# Assemble Figure 1d
g.fig1.1.corr <- (
  g.base[["corr"]][["sounds"]] + 
    (g.diff[["corr"]][["sounds"]] + plot_layout(tag_level = "new")) + 
    plot_layout(widths = c(2, 1), guides = "collect") &
    theme(strip.text = element_blank(),
          legend.position = "bottom",
          legend.margin = margin(t = -20, l = 60))
)

g.fig1.1.corr

Source code
# Assemble Supplementary Figure 1a
g.supfig2.corr <- (
  g.base[["corr"]][["omissions"]] + 
    (g.diff[["corr"]][["omissions"]] + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1)) &
    theme(strip.text.x = element_text(face = "bold"))
)

g.supfig2.corr

Theoretical null: accuracy plots

Source code
## Aggregate all data and clean

df.plot.maths <- bind_rows( 
  # --- make datasets compatible
  df.maths.mean %>% mutate(manip = "maths") %>% 
    mutate(t_train = dt*(t_train-1) - 0,
           t_test = dt*(t_test-1) - time_offset),
  df.stim.mean %>% 
    filter(manip == "") %>% 
    mutate(manip = "original") %>% 
    mutate(entropy = toupper(sub("rd_to_", "", direction))) %>% 
    # mutate(across(starts_with("t_"), ~(dt*(.-1)-700))) %>% 
    select(-stim, -direction)
) %>% 
  # --- calculate difference
  pivot_wider(names_from = manip, values_from = "accuracy") %>%
  mutate(diff = original - maths) %>%
  pivot_longer(cols = c(original, maths, diff),
               names_to = "manip", values_to = "accuracy") %>% 
  # --- crop time window
  filter(t_train>=0, (t_train-333) < dt) %>% 
  # --- aggregate results of cluster analysis
  left_join(df.clusters.maths) %>% 
  # --- clean labels
  mutate(manip = (recode_factor(manip, original = "REPLICATION\n", maths = "THEORETICAL NULL\n",
                                diff = "REPLICATION – THEORETICAL NULL\n")),
         direction = factor(entropy, levels = rev(c("RD","MM","MP","OR"))))
Joining with `by = join_by(t_train, t_test, entropy, manip)`
Source code
## PLOT
# Initialize subplots
g.maths.base <- list()
g.maths.diff <- list()

# --- ORIGINAL & THEORETICAL NULL
df.plot.maths %>% 
  filter(manip != "REPLICATION – THEORETICAL NULL\n") %>% 
  plot_base(col.pal = col.pal.base,
            z.breaks = acc.breaks.base,
            z.labels = acc.labels.base) +
  facet_grid(direction ~ manip) -> g

g.maths.base[["acc"]] <- g %>% add_onset(df.labels.onset[["sounds"]] %>% 
                                           mutate(manip = str_replace(manip, "EMPIRICAL", "THEORETICAL")))

# --- DIFFERENCE
df.plot.maths %>%
  filter(manip == "REPLICATION – THEORETICAL NULL\n") %>% 
  plot_base(col.pal = col.pal.diff,
            z.breaks = acc.breaks.diff,
            z.labels = acc.labels.diff) +
  facet_grid(direction ~ manip) +
  theme(axis.title.y = element_blank()) -> g

g.maths.diff[["acc"]] <- g %>% add_onset(df.labels.onset[["sounds"]] %>% 
                                           mutate(manip = str_replace(manip, "EMPIRICAL", "THEORETICAL")))

g.fig1.2.acc <- (
  g.maths.base[["acc"]] + 
    (g.maths.diff[["acc"]]  + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1)) &
    theme(strip.text.x = element_text(face = "bold"))
  )

g.fig1.2.acc

Theoretical null: correlations plots

Source code
df.plot <- df.plot.maths.corr %>% 
  mutate(manip = str_replace(manip, "ORIGINAL", "REPLICATION"),
         manip = paste0(manip, "\n"))

# ORIGINAL & THEORETICAL
df.plot %>% 
  filter(manip != "REPLICATION – THEORETICAL NULL\n") %>% 
  mutate(direction = "") %>% 
  plot_base(col.pal = col.pal.corr,
            z = "r",
            z.breaks = corr.breaks,
            z.labels = corr.labels) +
  labs(fill = "Correlation") +
  facet_grid(direction ~ manip) -> g

g.maths.base[["corr"]] <- g

# DIFFERENCE
df.plot %>%
  filter(manip == "REPLICATION – THEORETICAL NULL\n") %>% 
  mutate(direction = "") %>% 
  plot_base(col.pal = col.pal.corr,
            z = "r",
            z.breaks = corr.breaks,
            z.labels = corr.labels) +
  facet_grid(. ~ manip) +
  labs(fill = "Correlation") +
  theme(axis.title.y = element_blank()) -> g


g.maths.diff[["corr"]] <- g

g.fig1.2.corr <- (
  g.maths.base[["corr"]] + 
    (g.maths.diff[["corr"]] + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1), guides = "collect") &
    theme(strip.text.x = element_blank(),
          legend.position = "bottom",
          legend.margin = margin(t = -20, l = 60))
  )

g.fig1.2.corr

Assemble Figure 1

Source code
((g.transitions + theme(panel.spacing.x = unit(5,"mm"), legend.box.margin = margin(t = -30))) +
   g.reordering + theme(axis.text.x = element_text(margin = margin(t = -15))) +
   plot_layout(widths = c(1.7, 1))) / 
  ((g.fig1.1.acc / g.fig1.1.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) + 
     plot_layout(heights = c(4.7, 1))) / 
  ((g.fig1.2.acc / g.fig1.2.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) + 
     plot_layout(heights = c(4.7, 1))) +
  plot_layout(heights = c(1., 4, 4)) +
  plot_annotation(tag_levels = "a") -> g.fig1

g.fig1

Source code
if (save) {
  g.fig1 %>% ggsave(filename = "./figures/fig1_MattersArising.pdf", width = 7, height = 11, device = cairo_pdf)
}

Assemble Supplementary Figure 2 (omissions)

Source code
((g.supfig2.acc / g.supfig2.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) + 
   plot_layout(heights = c(4, 1))) +
  plot_annotation(tag_levels = "a") -> g.supfig2

g.supfig2

Source code
if (save) {
  g.supfig2 %>% ggsave(filename = "./figures/supfig2_MattersArising.pdf", width = 7, height = 5.5, device = cairo_pdf)
  g.supfig2 %>% ggsave(filename = "./figures/supfig2_MattersArising.png", width = 7, height = 5.5, device = ragg::agg_png(), dpi = 600)
}

Figure 2

Initialize

Source code
g.base.ml <- list()
g.diff.ml <- list()

Build accuracy plots

Source code
for (s in c("sounds")) {
  # ORIGINAL & REORDERED
  df.plot.ml %>% 
    filter(t_test >=-333, t_test <= 0) %>% 
    filter(stim == s, manip != "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip, 
                                 ORIGINAL = "ORIGINAL", 
                                 REORDERED = "EMPIRICAL NULL")) %>%
    plot_base(col.pal = col.pal.base,
              z.breaks = acc.breaks.base,
              z.labels = acc.labels.base,
              legend.position = "right") +
    facet_grid(direction ~ manip) +
    labs(title = paste0("Decoding MOST LIKELY sounds")) +
    theme(plot.title = element_text(hjust = -0.2)) -> g.base.ml[["acc"]][[s]]
  
  # DIFFERENCE
  df.plot.ml %>%
    filter(t_test >=-333, t_test <= 0) %>% 
    filter(!grepl("AVERAGE",direction)) %>% 
    filter(stim == s, manip == "ORIGINAL – REORDERED") %>% 
    mutate(manip = "DIFF") %>% 
    plot_base(col.pal = col.pal.diff,
              z.breaks = acc.breaks.diff,
              z.labels = acc.labels.diff,
              legend.position = "right") +
    facet_grid(direction ~ manip) +
    theme(axis.title.y = element_blank()) -> g.diff.ml[["acc"]][[s]]
}

(g.base.ml[["acc"]][["sounds"]] + (g.diff.ml[["acc"]][["sounds"]] + plot_layout(tag_level = "new")) + plot_layout(widths = c(2, 1))) + 
  theme(plot.tag.position = c(0, 1)) -> g.fig2.acc

g.fig2.acc
Warning: `stat_contour()`: Zero contours were generated
Warning in min(x): no non-missing arguments to min; returning Inf
Warning in max(x): no non-missing arguments to max; returning -Inf

Build correlations plots

Source code
for (s in c("sounds")) {
  # ORIGINAL & REORDERED
  df.plot.ml.corr %>% 
    filter(t_test >=-333, t_test <= 0) %>% 
    filter(stim == s, manip != "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip, 
                                 ORIGINAL = "ORIGINAL", 
                                 REORDERED = "EMPIRICAL NULL")) %>%
    mutate(direction = "") %>% 
    plot_base(col.pal = col.pal.corr,
              z = "r",
              z.breaks = corr.breaks,
              z.labels = corr.labels,
              legend.position = "right") +
    labs(fill = "Correlation") +
    facet_grid(direction ~ manip) -> g.base.ml[["corr"]][[s]]
  
  # DIFFERENCE
  df.plot.ml.corr %>%
    filter(t_test >=-333, t_test <= 0) %>% 
    filter(stim == s, manip == "ORIGINAL – REORDERED") %>%
    mutate(manip = "DIFF") %>% 
    mutate(direction = "") %>%
    plot_base(col.pal = col.pal.corr,
              z = "r",
              z.breaks = corr.breaks,
              z.labels = corr.labels,
              legend.position = "right") +
    facet_grid(. ~ manip) +
    labs(fill = "Correlation") +
    theme(axis.title.y = element_blank()) -> g.diff.ml[["corr"]][[s]]
}

(g.base.ml[["corr"]][["sounds"]] + g.diff.ml[["corr"]][["sounds"]] + plot_layout(widths = c(2, 1))) +
  # plot_annotation(tag_levels = list("b.")) &
  theme(plot.tag.position = c(0, 1)) -> g.fig2.corr

g.fig2.corr
Warning: `stat_contour()`: Zero contours were generated
Warning in min(x): no non-missing arguments to min; returning Inf
Warning in max(x): no non-missing arguments to max; returning -Inf
Warning: `stat_contour()`: Zero contours were generated
Warning in min(x): no non-missing arguments to min; returning Inf
Warning in max(x): no non-missing arguments to max; returning -Inf

Assemble figure 2

Source code
(g.mostlikely + coord_fixed(0.5))+
  plot_spacer() +
  ((g.fig2.acc / (g.fig2.corr + plot_layout(tag_level = 'new'))) + 
     plot_layout(heights = c(3, 1)) & 
     theme(plot.title = element_text(vjust = -8),
           legend.box.margin = margin(l=-5),
           strip.clip = "off",
           strip.text = element_text(size = 8),
           axis.text.x = element_text(angle = 45, hjust = 1))
  ) + 
  plot_layout(widths = c(1, 0.1, 2)) +
  plot_annotation(tag_levels = "a") & 
  theme(plot.tag.position = c(0, 1),
        plot.tag = element_text(vjust = 3)) -> g.fig2

g.fig2

Source code
if (save) {
  g.fig2 %>% ggsave(filename = "./figures/fig2_MattersArising.pdf", width = 8, height = 5.5, device = cairo_pdf)
}

Supplementary Figure 1

Prepare data

Source code
# Load subject-level data (used for clustering) and derive grand average (used for visualization) **without cropping**
tic()
df.supfig1.mean <- tibble()
df.supfig1.subj <- tibble()
for (manip in c("","_reord")) { 
  
  for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
    tmp <- load_loop_subj(folder = "sounds",
                          condition = paste0(direction, manip),
                          crop = FALSE)
    
    df.supfig1.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
    df.supfig1.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  }
}
toc() # ~1s
1.195 sec elapsed

Plot

Source code
# Extract the diagonal elements from the temporal generalization matrix
tmp <- filter(df.supfig1.subj, direction == "rd_to_rd", manip == "")[[1, "array"]][[1]]

tmp.diag <- matrix(nrow = dim(tmp)[3], ncol = dim(tmp)[1])
for (i in 1:dim(tmp)[3]) {
  tmp.diag[i,] <- tmp[,,i] %>% diag
}

# Plot
df.supfig1.mean %>% 
  filter(t_train == t_test,
         direction == "rd_to_rd",
         manip == "") %>% 
  mutate(sd = tmp.diag %>% apply(MARGIN = 2, FUN = sd) %>% `/`(sqrt(33))) %>% 
  filter(t_train >= -300, t_train <= 700) %>%
  
  ggplot(aes(x = t_train, y = accuracy)) +
  geom_hline(yintercept = 0.25, color = "grey50", linetype = 2) +
  geom_ribbon(aes(ymin = accuracy-sd, ymax = accuracy+sd), alpha = .3) +
  geom_line(size = 0.8) + 
  scale_x_continuous(breaks = scales::pretty_breaks(10), name = "Time (ms)") +
  scale_y_continuous(breaks = seq(.20, .40, .02), limits = c(.2,.4), name = "Accuracy") +
  coord_cartesian(expand = F)-> g.supfig1

g.supfig1

Source code
if (save) {
  g.supfig1 %>% ggsave(filename = "./figures/supfig1_MattersArising.pdf", width = 4, height = 2.5, device = cairo_pdf)
  g.supfig1 %>% ggsave(filename = "./figures/supfig1_MattersArising.png", width = 4, height = 2.5, device = ragg::agg_png(), dpi = 600)
}

Supplementary Figure 3

Functions & parameters

Source code
# Function to calculate and extract BF
getBF01 <- function(y, prior = "medium") {
  ttestBF(y, rscale = prior) %>% extractBF(onlybf = T)
}

# Color palette for BF
col.pal.bf <- rev(pals::brewer.piyg(7)); col.pal.bf[c(4)] <- "grey50";

Prepare data

Source code
tic()

## BF for EMPIRICAL null, DIFF
# --- extract array of subject-level correlation between accuracy and entropy 
X <- df.corr.subj[[3,"array"]][[1]]
# --- calculate BF for each cell of the TG matrix
bf.corr.diff.empir <- apply(X, c(1,2), getBF01, simplify = TRUE)


## BF for THEORETICAL null, DIFF
# --- extract array of subject-level correlation between accuracy and entropy 
X <- df.corr.maths.subj[[3,"array"]][[1]]
# --- calculate BF for each cell of the TG matrix
bf.corr.diff.maths <- apply(X, c(1,2), getBF01, simplify = TRUE)

## BF for MOST LIKELY null, DIFF
# --- extract array of subject-level correlation between accuracy and entropy 
X <- df.corr.ml.subj[[3,"array"]][[1]]
# --- calculate BF for each cell of the TG matrix
bf.corr.diff.ml <- apply(X, c(1,2), getBF01, simplify = TRUE)

toc() # ~1.5min
87.212 sec elapsed
Source code
## Post-processing
df.plot <- bind_rows(
  # --- aggregate data
  left_join(matrix_to_df2(bf.corr.diff.empir, "t_train", "t_test", "BF"), # BF results
            filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
    mutate(approach = "EMPIRICAL\nNULL"),
  left_join(matrix_to_df2(bf.corr.diff.maths, "t_train", "t_test", "BF"), # BF results
            filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
    mutate(approach = "THEORETICAL\nNULL"),
  left_join(matrix_to_df2(bf.corr.diff.ml, "t_train", "t_test", "BF"), # BF results
            filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
    mutate(approach = "MOST LIKELY")) %>% 
  
  # --- order factor levels 
  mutate(approach = factor(approach, levels = c("EMPIRICAL\nNULL", "THEORETICAL\nNULL", "MOST LIKELY"))) %>% 
           
  # --- convert indices to time
  mutate(t_test = (dt * (t_test - 1) - time_offset)) %>% 
  mutate(t_train = (dt * (t_train - 1))) %>% 
  # --- discretized version of BF
  mutate(BFd = cut(BF,
                   breaks = c(0, 1/30, 1/10, 1/3, 3, 10, 30, Inf),
                   labels = c("very strong in favor of null (BF < 1/30)",
                              "strong in favor of null (1/30 < BF < 1/10)",
                              "moderate in favor of null (1/10 < BF < 1/3)",
                              "inconclusive (1/3 < BF < 3)",
                              "moderate in favor of alt. (3 < BF < 10)",
                              "strong in favor of alt. (10 < BF < 30)",
                              "very strong in favor of alt. (BF > 30)")),
         .after = 3
         ) %>% 
  # --- transform continuous BF for easy symetric display on colorbar
  mutate(BF = case_when(BF < 1/3 ~ -1/BF,
                        between(BF, 1/3, 3) ~ NA,
                        BF > 3 ~ BF))
Joining with `by = join_by(t_train, t_test)`
Joining with `by = join_by(t_train, t_test)`
Joining with `by = join_by(t_train, t_test)`
Source code
df.plot %>% summary
    t_train        t_test           BF            
 Min.   :  0   Min.   :-670   Min.   :-5.000e+00  
 1st Qu.: 80   1st Qu.:-340   1st Qu.:-5.000e+00  
 Median :160   Median :   0   Median :-5.000e+00  
 Mean   :160   Mean   :   0   Mean   : 2.063e+06  
 3rd Qu.:240   3rd Qu.: 340   3rd Qu.:-4.000e+00  
 Max.   :320   Max.   : 670   Max.   : 1.539e+10  
                              NA's   :4029        
                                          BFd           signif      
 very strong in favor of null (BF < 1/30)   :   0   Min.   :0.0000  
 strong in favor of null (1/30 < BF < 1/10) :   0   1st Qu.:0.0000  
 moderate in favor of null (1/10 < BF < 1/3):7945   Median :0.0000  
 inconclusive (1/3 < BF < 3)                :4029   Mean   :0.2496  
 moderate in favor of alt. (3 < BF < 10)    : 569   3rd Qu.:0.0000  
 strong in favor of alt. (10 < BF < 30)     : 264   Max.   :1.0000  
 very strong in favor of alt. (BF > 30)     : 558                   
    manip               stim                      approach     
 Length:13365       Length:13365       EMPIRICAL\nNULL  :4455  
 Class :character   Class :character   THEORETICAL\nNULL:4455  
 Mode  :character   Mode  :character   MOST LIKELY      :4455  
                                                               
                                                               
                                                               
                                                               

Plot

Source code
lines.v.pos <- 330*seq(round(min(df.plot$t_test)/333), round(max(df.plot$t_test)/333), 1)
lines.v.pos <- lines.v.pos[lines.v.pos!=0]
lines.h.pos <- 330*seq(round(min(df.plot$t_train)/333), round(max(df.plot$t_train)/333), 1)
lines.h.pos <- lines.h.pos[lines.h.pos!=0]


df.plot %>%
  ggplot(aes(x = t_test, y = t_train)) +
  facet_grid(approach ~ .) +
  geom_tile(aes(fill = BFd), show.legend = TRUE) +
  
  geom_vline(xintercept = lines.v.pos, linetype = 2, color = "black") +
  geom_hline(yintercept = lines.h.pos, linetype = 2, color = "black") +
  geom_vline(xintercept = 0, color = "black") +
  geom_hline(yintercept = 0, color = "black") +
  
  geom_contour(aes(z = signif), size = 0.2, color = col.clusters) +
  
  scale_x_continuous(breaks = scales::pretty_breaks(12),
                     expand = c(0,0)) +
  scale_y_continuous(breaks = scales::pretty_breaks(4),
                     expand = c(0,0)) +
  scale_fill_manual(values = rev(col.pal.bf),
                    labels = levels(df.plot$BFd),
                    drop = FALSE) +
  
  labs(fill = expression(BF[10])) +
  coord_equal() +
  
  labs(x = "Test time (ms)", y = "Train time (ms)") -> g.supfig3

g.supfig3

Source code
if (save) {
  g.supfig3 %>% ggsave(filename = "./figures/supfig3_MattersArising.pdf", width = 10, height = 3, device = cairo_pdf)
  g.supfig3 %>% ggsave(filename = "./figures/supfig3_MattersArising.png", width = 10, height = 3, device = ragg::agg_png(), dpi = 600)
}

Supplementary Figure 4

Source code
mat.C.all <- array(dim = c(34,34,4,4,n.subj))
i <- 1

for (s in list.subj) {
    
    # --- load the confusion matrix
    tmp.file <- os$path$join(path.root, s, "rd_to_rd_confmats.npz") %>%
      np$load()
    
    # --- crop around S0 (decoded stimulus)
    mat.C.all[,,,,i] <- tmp.file$f[["arr_0"]][,71:104,71:104,,] %>%
      apply(c(2,3,4,5), mean)
    
    i <- i+1
}

mat.C.mean <- mat.C.all %>% apply(c(1,2,3,4), mean)

    
df.confusion <- matrix_to_df2(mat.C.mean[22,22,,], "from", "to", "p") %>% 
  mutate(label = paste0("a[",from,to,"]"))
df.confusion %>% 
  ggplot(aes(x = to, y = from, fill = p)) +
  geom_tile(color = "firebrick3") + 
  geom_text(aes(label = label), parse = T, color = "firebrick3") +
  scale_x_continuous(name = "predicted class", position = "top") +
  scale_y_reverse(name = "true class") +
  scale_fill_stepsn(name = "p",
                    colors = pals::brewer.greys(100),
                      # low = "white", high = "black",
                      labels = ~str_replace(., "^0.","."),
                      breaks = scales::pretty_breaks(6),
                      limits = c(.20,.30), oob = scales::squish) +
  coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
              expand = F, clip = "off") -> g.matC

matrix_to_df2(mat.T$MP, "from", "to", "p") %>% 
  mutate(label = paste0("a[",from,to,"]")) %>% 
  ggplot(aes(x = to, y = from, fill = p)) +
  geom_tile(color = "firebrick3") + 
  geom_text(aes(label = label), parse = T, color = "firebrick3") +
  scale_x_continuous(name = "to", position = "top") +
  scale_y_reverse(name = "from") +
  scale_fill_stepsn(name = "p",
                    colors = pals::brewer.greys(100),
                      labels = ~str_replace(., "^0.","."),
                      breaks = scales::pretty_breaks(6),
                      limits = c(.0,.50), oob = scales::squish) +
  coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
              expand = F, clip = "off") -> g.matT

 
((g.matT + labs(caption = "Transition matrix\n(midplus (MP) sequence)") +
    guides(fill = guide_colorbar(barwidth = 0.5, ticks = F))) +
(g.matC + labs(caption = "Confusion matrix\n(group average, t=210ms)") +
    guides(fill = guide_colorbar(barwidth = 0.5, ticks = F))) &
  theme(plot.caption = element_text(size = 10),
        axis.title = element_text(face = "italic"),
        legend.box.spacing = unit(2,"mm"),
        )
  ) +
  plot_annotation(tag_levels = "a") -> g.supfig4

g.supfig4

Source code
if (save) {
  g.supfig4 %>% ggsave(filename = "./figures/supfig4_MattersArising.pdf", width = 7, height = 2.5, device = cairo_pdf)
  g.supfig4 %>% ggsave(filename = "./figures/supfig4_MattersArising.png", width = 7, height = 2.5, device = ragg::agg_png(), dpi = 600)
}

References

Demarchi, Gianpaolo, Gaëtan Sanchez, and Nathan Weisz. 2019. “Automatic and Feature-Specific Prediction-Related Neural Activity in the Human Auditory System.” Nature Communications 10 (1): 1–11. https://doi.org/10.1038/s41467-019-11440-1.